# 计算特征空间
import numpy as np
from matplotlib import pyplot as plt


def feature_space(X, Y):
    """
    :param X: 自变量
    :param Y: 因变量
    :return: 特征空间拟合系数
    """
    # 对输入数据进行异常值剔除，最好在输入前剔除，这里仅做通用剔除
    X[np.isinf(X)] = np.nan
    Y[np.isinf(Y)] = np.nan
    # 归一化
    # X = (X - np.nanmean(X)) / (np.nanmax(X) - np.nanmin(X))
    # Y = (Y - np.nanmean(Y)) / (np.nanmax(Y) - np.nanmin(Y))
    # 第一步，将evi划分为M个区间（建议M>=20)，每个区间分成N个子区间（建议N>=5)
    M = 20
    N = 5
    # 制定自变量取值一级区间
    x_interval = np.linspace(np.nanmin(X), np.nanmax(X), M + 1, endpoint=True)
    print(x_interval)
    # 第二步，对给定区间找出每个子区间的最高因变量和最低因变量保存
    y_max_list = []
    y_min_list = []
    # 满足条件的自变量节点
    x_max_list = []
    x_min_list = []
    for i in range(M):
        # 制定子区间取值
        x_sub_interval = np.linspace(x_interval[i], x_interval[i + 1], N + 1, endpoint=True)
        # 存储子区间中的最大值和最小值
        y_sub_max_list = []
        y_sub_min_list = []
        for j in range(N):
            x_sub_interval_left = x_sub_interval[j]
            x_sub_interval_right = x_sub_interval[j + 1]
            # 获取在这个子区间里面的因变量值
            y_sub_interval_value = Y[np.where((X >= x_sub_interval_left) & (X < x_sub_interval_right))]
            # print(y_sub_interval_value)
            if y_sub_interval_value.size != 0:
                y_sub_max = np.nanmax(y_sub_interval_value)
                y_sub_min = np.nanmin(y_sub_interval_value)
                if ~np.isnan(y_sub_max):
                    y_sub_max_list.append(y_sub_max)
                if ~np.isnan(y_sub_min):
                    y_sub_min_list.append(y_sub_min)
        # 判断该区间是否有值，没有值则跳过
        if y_sub_max_list:
            need_remove = True
            y_max_mean = np.nan
            # 第三步，计算给定区间的N个子区间的标准偏差，直到满足偏差
            while need_remove:
                y_max_mean = np.nanmean(y_sub_max_list)
                # 如果里面不为nan的元素少于两个就不按照偏差剔除了
                if len(y_sub_max_list) > 2:
                    for j in range(len(y_sub_max_list)):
                        y_max_std = np.nanstd(y_sub_max_list)
                        # 第四步，如果给定区间的每个子区间的最高温度小于平均值-偏差，则在下面步骤丢弃该子区间
                        # todo 待商榷
                        if y_sub_max_list[j] < (y_max_mean - y_max_std) and ~np.isnan(y_sub_max_list[j]):
                            y_sub_max_list.remove(y_sub_max_list[j])
                            need_remove = True
                            break
                        else:
                            need_remove = False
                else:
                    break
            if ~np.isnan(y_max_mean):
                y_max_list.append(y_max_mean)
                x_max_list.append(x_interval[i])
        # print("y_sub_min_list: ", y_sub_min_list)
        if y_sub_min_list:
            need_remove = True
            y_min_mean = np.nan
            # 第三步，计算给定区间的N个子区间的标准偏差，直到满足偏差
            while need_remove:
                y_min_mean = np.nanmean(y_sub_min_list)
                # 如果里面元素少于两个就不按照偏差剔除了
                if len(y_sub_min_list) > 2:
                    for j in range(len(y_sub_min_list)):
                        y_min_std = np.nanstd(y_sub_min_list)
                        # 第四步，如果给定区间的每个子区间的最高温度小于平均值-偏差，则在下面步骤丢弃该子区间
                        # todo 待商榷
                        if y_sub_min_list[j] < (y_min_mean - y_min_std) and ~np.isnan(y_sub_min_list[j]):
                            y_sub_min_list.remove(y_sub_min_list[j])
                            need_remove = True
                            break
                        else:
                            need_remove = False
                else:
                    break
            if ~np.isnan(y_min_mean):
                y_min_list.append(y_min_mean)
                x_min_list.append(x_interval[i])
    print(x_min_list)
    print(y_min_list)
    print(x_max_list)
    print(y_max_list)
    # 第六步，如果给定区间内剩余子区间数量大于给定阈值且标准差大于给定阈值，则返回步骤四，重复步骤四-六。否则转到步骤七
    # 看是否给阈值，这里暂时给了2个元素
    # 第七步，取平均温度作为该区间的最高温度，并返回步骤二，知道找到M个区间内所有最高温度。
    # 前面循环已做
    # 第八步，线性回归，并计算RMSE。
    a, b, c, d = 0, 0, 0, 0
    # 拟合干边
    need_fit = True
    while need_fit:
        y_max_arr = np.array(y_max_list)
        x_max_arr = np.array(x_max_list)
        popt = np.polyfit(x_max_arr, y_max_arr, 1)
        a = popt[0]
        b = popt[1]
        # print(a, b)
        y_max_fit = a * x_max_arr + b
        y_max_fit_rmse = np.sqrt(np.mean((y_max_arr - y_max_fit) ** 2))
        if len(y_max_list) > 2:
            for i in range(len(y_max_list)):
                # 第九步，如果给定区间最高温度是RMSE的2倍或大于等于回归线上的温度值，则该区间将被丢失，返回步骤八，直到达到最小区间数或不能在丢弃任何区间为止
                if y_max_list[i] < y_max_fit[i] - y_max_fit_rmse * 2:
                    y_max_list.remove(y_max_list[i])
                    need_fit = True
                    break
                else:
                    need_fit = False
        else:
            break

    # 拟合湿边
    need_fit = True
    while need_fit:
        y_min_arr = np.array(y_min_list)
        x_min_arr = np.array(x_min_list)
        popt = np.polyfit(x_min_arr, y_min_arr, 1)
        c = popt[0]
        d = popt[1]
        # print(c, d)
        y_min_fit = c * x_min_arr + d
        y_min_fit_rmse = np.sqrt(np.mean((y_min_arr - y_min_fit) ** 2))
        if len(y_min_list) > 2:
            for i in range(len(y_min_list)):
                # 第九步，如果给定区间最高温度是RMSE的2倍或大于等于回归线上的温度值，则该区间将被丢失，返回步骤八，直到达到最小区间数或不能在丢弃任何区间为止
                if y_min_list[i] < y_min_fit[i] - y_min_fit_rmse * 2:
                    y_min_list.remove(y_min_list[i])
                    need_fit = True
                    break
                else:
                    need_fit = False
        else:
            break

    # 第十步，执行最后的线性回归得到干边，
    y_max_interval = a * x_interval + b
    y_min_interval = c * x_interval + d
    plt.title("NDVI-ATI")
    plt.scatter(X.flatten(), Y.flatten(), c='r', marker='o')
    plt.plot(x_interval, y_max_interval, c='b', label='max')
    plt.plot(x_interval, y_min_interval, c='g', label='min')
    plt.xlabel("NDVI")
    plt.ylabel("ATI")
    plt.legend()
    plt.show()
    return a, b, c, d